# Script 1 - raw vulnerability curve data
# Figure 1

#################################################
##### Script to analyse raw optical data ########
#################################################
# RP Skelton 
# July 2019
# edited Nov 2022
setwd("C://Users/rober/Dropbox/2024_ThicketHydraulics_SkeltonButtnerPotts/1_data")

setwd("D:/Dropbox/100_PROJECTS/2024_ThicketHydraulics_SkeltonButtnerPotts/1_data")

# load functions
se <- function(x) sqrt(var(x, na.rm=T)/sum(!is.na(x)))

# load optical data
all.ovc <- read.csv("04_OpticalVulnerabilityCurves.csv",header=T)
l <- 0
# Create ov.data dataframe
ov.data <- matrix(NA,nrow = 30,ncol=7)
ov.data <- as.data.frame(ov.data)
colnames(ov.data) <- c("genus","species","individual","pe","p50","p12","slope")

# find summary data and add to figure

##############################################
########## Generating VC data ################
##############################################
cols <- c("#F2EDD8","#F2EDD8","#F2EDD8")

# Schotia latifolia
# Individual 1
x <- all.ovc$WP[all.ovc$Genus  == "Schotia" & all.ovc$Individual == 1]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 1]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Schotia"
ov.data[l,2] <- "latifolia"
ov.data[l,3] <- 1
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
ov.data[l,7] <- coefficients(summary(plcmod))[1]
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 2
x <- all.ovc$WP[all.ovc$Genus  == "Schotia" & all.ovc$Individual == 2]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 2]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Schotia"
ov.data[l,2] <- "latifolia"
ov.data[l,3] <- 2
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 3
x <- all.ovc$WP[all.ovc$Genus  == "Schotia" & all.ovc$Individual == 4]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Schotia" & all.ovc$Individual == 4]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Schotia"
ov.data[l,2] <- "latifolia"
ov.data[l,3] <- 3
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 4
x <- all.ovc$WP[all.ovc$Genus  == "Schotia" & all.ovc$Individual == 5]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 5]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Schotia"
ov.data[l,2] <- "latifolia"
ov.data[l,3] <- 4
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]

# Polygala myrtifolia
# Individual 1
x <- all.ovc$WP[all.ovc$Genus  == "Polygala" & all.ovc$Individual == 1]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 1]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Polygala"
ov.data[l,2] <- "myrtifolia"
ov.data[l,3] <- 1
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]

#Individual 2
x <- all.ovc$WP[all.ovc$Genus  == "Polygala" & all.ovc$Individual == 2]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 2]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Polygala"
ov.data[l,2] <- "myrtifolia"
ov.data[l,3] <- 2
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 3
x <- all.ovc$WP[all.ovc$Genus  == "Polygala" & all.ovc$Individual == 3]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 3]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Polygala"
ov.data[l,2] <- "myrtifolia"
ov.data[l,3] <- 3
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 4
x <- all.ovc$WP[all.ovc$Genus  == "Polygala" & all.ovc$Individual == 4]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 4]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Polygala"
ov.data[l,2] <- "myrtifolia"
ov.data[l,3] <- 4
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]

# Pappea capensis
# Individual 5
x <- all.ovc$WP[all.ovc$Genus  == "Pappea" & all.ovc$Individual == 5]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Pappea"& all.ovc$Individual == 5]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-5))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Pappea"
ov.data[l,2] <- "capensis"
ov.data[l,3] <- 1
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 2
x <- all.ovc$WP[all.ovc$Genus  == "Pappea" & all.ovc$Individual == 2]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Pappea"& all.ovc$Individual == 2]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-5))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Pappea"
ov.data[l,2] <- "capensis"
ov.data[l,3] <- 2
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 3
x <- all.ovc$WP[all.ovc$Genus  == "Pappea" & all.ovc$Individual == 3]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Pappea"& all.ovc$Individual == 3]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-5))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Pappea"
ov.data[l,2] <- "capensis"
ov.data[l,3] <- 3
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2]
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 4
x <- all.ovc$WP[all.ovc$Genus  == "Pappea" & all.ovc$Individual == 4]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Pappea"& all.ovc$Individual == 4]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Pappea"
ov.data[l,2] <- "capensis"
ov.data[l,3] <- 4
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]

# Boscia oleoides
# Individual 1
x <- all.ovc$WP[all.ovc$Genus  == "Boscia" & all.ovc$Individual == 1]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 1]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Boscia"
ov.data[l,2] <- "oleoides"
ov.data[l,3] <- 1
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]

#Individual 2
x <- all.ovc$WP[all.ovc$Genus  == "Boscia" & all.ovc$Individual == 2]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 2]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Boscia"
ov.data[l,2] <- "oleoides"
ov.data[l,3] <- 2
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 3
x <- all.ovc$WP[all.ovc$Genus  == "Boscia" & all.ovc$Individual == 3]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 3]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Boscia"
ov.data[l,2] <- "oleoides"
ov.data[l,3] <- 3
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 4
x <- all.ovc$WP[all.ovc$Genus  == "Boscia" & all.ovc$Individual == 4]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 4]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-4))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Boscia"
ov.data[l,2] <- "oleoides"
ov.data[l,3] <- 4
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]

# Euclea undulata
# Individual 1
x <- all.ovc$WP[all.ovc$Genus  == "Euclea" & all.ovc$Individual == 1]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Euclea"& all.ovc$Individual == 1]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-6))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Euclea"
ov.data[l,2] <- "undulata"
ov.data[l,3] <- 1
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 2
x <- all.ovc$WP[all.ovc$Genus  == "Euclea" & all.ovc$Individual == 2]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Euclea"& all.ovc$Individual == 2]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-6))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Euclea"
ov.data[l,2] <- "undulata"
ov.data[l,3] <- 2
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 3
x <- all.ovc$WP[all.ovc$Genus  == "Euclea" & all.ovc$Individual == 3]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Euclea"& all.ovc$Individual == 3]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-6))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Euclea"
ov.data[l,2] <- "undulata"
ov.data[l,3] <- 3
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2]
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]

# Searsia longispina
# Individual 1
x <- all.ovc$WP[all.ovc$Genus  == "Searsia" & all.ovc$Individual == 1]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Searsia"& all.ovc$Individual == 1]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-7))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Searsia"
ov.data[l,2] <- "longispina"
ov.data[l,3] <- 1
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 2
x <- all.ovc$WP[all.ovc$Genus  == "Searsia" & all.ovc$Individual == 2]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Searsia"& all.ovc$Individual == 2]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-7))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Searsia"
ov.data[l,2] <- "longispina"
ov.data[l,3] <- 2
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]
#Individual 3
x <- all.ovc$WP[all.ovc$Genus  == "Searsia" & all.ovc$Individual == 3]
y <- all.ovc$PercentEmbolism[all.ovc$Genus  == "Searsia"& all.ovc$Individual == 3]
plcmod <- nls(y~100-100/(1+exp(a*(x-b))), start=list(a=-6,b=-7))
xv.2 <- seq(from=-10,to=0,by=0.01)
yv.2 <- predict(plcmod,data.frame(x=xv.2))
l <- l + 1
ov.data[l,1] <- "Searsia"
ov.data[l,2] <- "longispina"
ov.data[l,3] <- 3
ov.data[l,4] <- mean(x[which(y == min(y[y>12]))],na.rm=T)
ov.data[l,5] <- coefficients(summary(plcmod))[2] 
### P12
a <- coefficients(summary(plcmod))[1]
p50 <- coefficients(summary(plcmod))[2]
slope <- -100*a/4
c <- (100/2) - slope*p50
p88 <- -c/slope
p12 <- (100-c)/slope
ov.data[l,6] <- p12
ov.data[l,7] <- coefficients(summary(plcmod))[1]

#########################################################
############### END OF ANALYSIS ONE #####################
############# Raw vulnerability curves ##################
#########################################################

#####################################################
##### ANOVA testing differences between species #####
#####################################################
# Conduct statistical tests to compare safety margins between species
# Tabulate the results
# P50
anova(lm(ov.data$p50~ov.data$species))
x <- aov(lm(ov.data$p50~ov.data$species))
summary(x)
TukeyHSD(x)
# boxplot(ov.data$p50~ov.data$species)
# Result:
# Some variation in P50 between species
# a  Schotia
# ab  Pappea
# a  Polygala
# b  Boscia
# c Searsia, Euclea

# Slope
anova(lm(ov.data$slope~ov.data$species))
x <- aov(lm(ov.data$slope~ov.data$species))
summary(x)
TukeyHSD(x)
# boxplot(ov.data$slope~ov.data$species)
# Result:
# Some variation in slope between species
# ab  Schotia
# a  Pappea
# a  Polygala
# a  Boscia
# a Euclea
# b Searsia

p50.mean <- round(tapply(ov.data$p50,ov.data$species,mean,na.rm=T),2)
p50.mean.se <- round(tapply(ov.data$p50,ov.data$species,se),2)

pe.mean.stem <- round(tapply(ov.data$pe,ov.data$species,mean,na.rm=T),2)
pe.mean.se.stem <- round(tapply(ov.data$pe,ov.data$species,se),2)

p12.mean.stem <- round(tapply(ov.data$p12,ov.data$species,mean,na.rm=T),2)
p12.mean.se.stem <- round(tapply(ov.data$p12,ov.data$species,se),2)

# slope
slope.mean <- round(tapply(ov.data$slope,ov.data$species,mean,na.rm=T),2)
slope.mean.se <- round(tapply(ov.data$slope,ov.data$species,se),2)

# sample sizes
ind.num <- tapply(ov.data$individual,ov.data$species,length)

summary.data <- matrix(NA, ncol=11, nrow=length(row.names(p50.mean)))
summary.data <- as.data.frame(summary.data)
colnames(summary.data) <- c("tip.label", "species", "P50", "P50se", "slope","slope.se","ss","pe","pe.se","p12","p12.se")
summary.data[,1] <- c("Pappea_capensis","Schotia_latifolia","Searsia_longispina","Polygala_myrtifolia","Boscia_oleoides","Euclea_undulata")
summary.data[,2] <- row.names(p50.mean)
summary.data[,3] <- p50.mean
summary.data[,4] <- p50.mean.se
summary.data[,5] <- slope.mean
summary.data[,6] <- slope.mean.se
summary.data[,7] <- (ind.num)[match(names(ind.num), summary.data$species)]
summary.data[,8] <- pe.mean.stem
summary.data[,9] <- pe.mean.se.stem
summary.data[,10] <- p12.mean.stem
summary.data[,11] <- p12.mean.se.stem
summary.data


write.csv(summary.data,"1_data/04b_OpticalVulnerabilityCurvesSummaryData.csv")